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How strong are the Rossby vortices? 
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ABSTRACT 

The Rossby wave instability, associated with density bumps in differentially rotating 
discs, may arise in several different astrophysical contexts, such as galactic or proto- 
planetary discs. While the linear phase of the instability has been well studied, the 
nonlinear evolution and especially the saturation phase remain poorly understood. In 
this paper, we test the non-linear saturation mechanism analogous to that derived for 
wave-particle interaction in plasma physics. To this end we perform global numerical 
simulations of the evolution of the instability in a two-dimensional disc. We confirm 
the physical mechanism for the instability saturation and show that the maximum 
amplitude of vorticity can be estimated as twice the linear growth rate of the instabil- 
ity. We provide an empirical fitting formula for this growth rate for various parameters 
of the density bump. We also investigate the effects of the azimuthal mode number 
of the instability and the energy leakage in the spiral density waves. Finally, we show 
that our results can be extrapolated to 3D discs. 
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1 INTRODUCTION 

Understanding the evolution of accretion discs is impor- 
tant for several astrop hysical applications. The Rossby 
wave instability (RWI) l| Lovelace et a.1.1 Il999l ) and Rossby 
vortices it form have been discussed in a large va- 
riety of disc systems, from large-scale galact ic discs 
l|Lovelace fc Hohlfeldll 19781; ISellwood fc Kahrjl99lf ) and the 
Galactic centre dTagger fc Melial l2006l h to discs in mi- 



croquasars jLovelace et al 



i|Varniere fc Tagged l200rl 



|2009|1 and protopla netary discs 



Meheut et al]|2012h . The latter 



are of particular importance since Rossby vortices may ac- 
celerate planetesimal formation. The RWI operates around 
the density bumps or vortensity bumps (see Section 2) in the 
disc. In the linear theory, the growth rate of the instability 
ca n be computed as an eigenvalue problem. This w as done 
bv lLi et alJ (|2000t l for 2D discs, bv lUmurhanl J20ld) in the 
shal low water approximation, and by I Meheut et all l|2012h 
and|Ei3 (|2012h for 3D vertically stratified discs. After the 
linear phase, the instability saturates and the linear theory 
breaks down. The saturation amplitude of the instability is 
of importance in the several applications. For example, in 
protoplanetary discs, the vortices formed by the RWI can 
concentrate the solid particles and hence accelerate the for- 
mation of planetesimals. The amplitude of the vortices is 
a key parameter that determines the amount of solids that 
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can be concentrated and t he size of the resulting planetes- 
imals (|Meheut et al.l [20121 ) . The Rossby vortices can also 
influence planet migration in the dis c in the latter stage of 
planet formation (|Koller et al.ll2003T ). and the amplitude of 
migration is determined by the strength of the vortices. In 
the RWI model of Quasi-Periodic Oscillations (QPO) in mi- 
croquasars, the saturation amplitude of the RWI is of impor- 
tance for the understanding of the observed QPO amplitude. 

The saturation mec hanism of the RWI was already dis- 
cussed in iLi et al.li|200llh but no criterion for saturation was 
given. In this paper, we interpret the growth of the insta- 
bility in term of particle-wave interaction and use classical 
plasma physics results in this context, whic h give the non- 
linear saturation condition as proposed by lLovelace et al.l 
(2009). In section 2, we present the physical interpretation 
of the growth of the instability and the model for the satu- 
ration mechanism. The numerical simulations are presented 
and compared with the model in section [3] We then sum- 
marise and conclude in the last section. 



2 THE RWI 

2.1 Linear growth 

The RWI can be seen as an equivalent of the Kelvin- 
Helmholtz instability in differentially rotating discs, and has 
a similar instability criterion: an extremum in the quantity 
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Figure 1. Schematic view of the RWI with the two propagating 
regions for the Rossby waves and density waves, and in between 
the evanescent regions. At the corotation radius r c the wave pat- 
tern frequency matches the orbital frequency (Eq. [4) . The inner 
and outer Lindblad radii, ?*ilr and rpLRi a- re defined by Eq. J5} 



£ related to the vorticity of the equilibrium flow. In a non- 
magnetised thin disc, this quantity can be written as 



(1) 



£ = ^(pE-^ = —J^ (pE-T) 

^ 2(V x v) z 

where E is the surface density, v the velocity of the fluid, 
7 the adiabatic index, Q. the rotation frequency, and k 2 = 
40 2 -\-2rQQ' the squared epicyclic frequency (so that « 2 /2f2 
is the vorticity). Here the prime denotes a radial derivative. 
For barotropic discs considered in this paper, the quantity 
£ reduces to 



C 



Efi 
~7T 



(2) 



2(V x v) z 

corresponding to 1/2 of the inverse vortensitjQ. 

In a disc with such an extremum, two standing Rossby 
waves can coexist, one on each side of the extremum. We re- 
call here that Rossby waves are vorticity waves propagating 
on the gradient of vortensity, and their dispersion relation 
is given by 



2k 2 ci 



mC 



c 2 {k 2 +m 2 /r 2 ) + K 2 rE 



(3) 



where u is the wave frequency, cD = cu — mQ is the wave 
frequency in the rotating frame, m and k r its azimuthal and 
radial wave numbers, and c s is the sound speed. The wave 
corotation radius r c is defined by: 



w(r c ) 



mQ,(r c 



0. 



(4) 



Thus, with the negative gradient of the disc rotation velocity, 
we have Cj r>rc < and Cj r <r c > 0. From the dispersion 
relation, one can see that Rossby waves propagate only in 
the region where £' > if r < r c and where £' < if r > r c . 

We consider here the case of a maximum of £ as plot- 
ted in Fig. [T] The wave located in the positive gradient of 
£ (r < r c ) has the pattern frequency smaller than the disc 
rotation frequency (co/m < SI) and carries negative energy 



Vortensity is defined as the ratio of vertical vorticity to surface 
density: ( v 



(i.e., increasing wave amplitude tends to to remove energy 
from the fluid) . On the other hand, the wave in the negative 
gradient of £ has uj/m > Q, and carries positive energy. The 
growth of the RWI is due to the interaction between the two 
Rossby waves with respectively positive and negative ener- 
gies. The total energy is conserved while the wave amplitude 
increases. The mechanism for the growth of the instability is 
then localised in the corotation region and does not depend 

on the boundary con ditions. 

As explained by Tagger (2001), in accretion discs the 
differential rotation couples Rossby waves to density waves. 
These density waves can only propagates outside the two 
(inner and outer) Lindblad resonances radii tlr defined as 



lo - mfi(rLR.) = ±«(rLR). 



(5) 



Between the two propagation zones the wave is evanescent 
and tunnelling is po ssible as plotted in Fig. [T] As detailed in 
iTsang fc LaH IMPOST ) and lLai fc Tsand (120091 ') . these propa- 
gation regions clearly appear when expressing the linearized 
perturbation equations in the form of Schrodinger equation, 
with the effective potential given by 



2mSl d 
r£j d r 



( ln ^^) 



+ 



+ 



(6) 



The waves are propagating in the region where V a s (r) is neg- 
ative and evanescent where V e s(r) is positive. This effective 
potential is plotted in Fig. [2] for the specific configuration 
considered below. 



2.2 Non-linear saturation 

To estimate the amplitude at which the linear calcula- 
tion stops to be valid, we use a standard plasma crite- 
rion which gives, for instance, the breakdow n of the Lan- 
dau damping theory due to particle trapping (|Q'Neillll965l : 
iKrall fc Trivelpiecdll9"73h . We consider a fluid particle in a 
Rossby wave with the vorticity amplitude uj v exponentially 
growing with time during the linear stage, 

\uh\ oc exp(Tf). (7) 

We here use the convention to express time in units of 

n -\ 

The estimation of the rotation time of a fluid particule 
in a vortex is related to the vorticity. For instance for a 
circular vortex with constant vorticity, it can be estimated 
as 



tt ~ 47r/|ai„ | 



(8) 



However in the case of the vortices formed by the RWI, 
the vortex doesn't have an analytical expression and the 
turnover timescale can not be computed analytically. From 
the simulations presented in the next section, we estimated 



tt 



2/\oj v 



(9) 



This circular motion around the vortex centre is a major 
perturbation to the circular orbit of the fluid particle around 
the disc centre and it modifies the radial structure of the disc 
in the region of the vortices on a time scale tt- The linear 
calculation of the growth of the Rossby wave instability is 
determined by this radial structure of the disc and especially 
by the radial structure of £. Therefore, this linear approach 
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is only valid when the growth timescale (I/7) is smaller that 
tt- This occurs when 



\L0 V 

~2~ 



7. 



(10) 



In other words, the exponential growth of the wave am- 
plitude necessarily ends when the fluid particle motion in 
the vicinity of resonant radius r c differs appreciably from 
the strictly azimuthal motion in the equilibrium. Note that 
the saturation criterion (|10p is onl y order-of-magn it ude, and 
similar criterion was proposed by lLovelace et al.l (j2009h in 
a different context. Below we perform non-linear numerical 
simulations to calibrate this criterion. 



3 NON-LINEAR SIMULATIONS 
3.1 Methods 

We have performed two dimensional (2D) non-linear simula- 
tions of a differentially rotating disc with a constant surface 
density except for a Gaussian bump. This over-density gives 
an extremum of C in which the instability can grow. 

We work in cylindrical coordinates (r, <p) with the 2D 
Euler equations 



<9tE + V-(ptT) = 0, 
di(EtT) + V ■ (vEv) + Vp = -EV$ G . 



(11) 
(12) 



where E is the surface density of the fluid, v its velocity, and 
p its pressure. <E>g = —GM,/r is the gravitational potential 
of the central object with G the gravitational constant and 
M» the mass of the star. We consider a globally isothermal 
flow, i.e. with a radially constant sound speed, and a linear 
relation between pressure and density p = c^E. 

The initial surface density, normalized to Eq is given by 



E/E = l + xexp(-^). 



(13) 



Our canonical parameters are r = 1., the amplitude of the 
bump x is chosen in the range [0.15,0.3], and its width is 
given by afro = 0.05. Although we will also present results 
for different values of a/ro. The disc scale height is fixed to 
h — c 3 /Qq = O.lro, S7q being the Keplerian orbital frequency 
at ro. The initial azimuthal velocity is chosen to achieve 
radial equilibrium. A low amplitude perturbation is added 
to this equilibrium: 



v T = esm(m<f)) exp I — 



(r - r a f 
2a 2 



(14) 



with e = 10~ 4 and the azimuthal mode number m £ [2,5]. 
We have also performed simulations with initial random per- 
turbations of the same amplitude. 

We use the Message Passing Interface-Adaptive Mesh 
Refinement Versatile Advection Code (MPI-AMRVAC) 
jKeppens et al.ll2012f ), with a uniform grid. The numeri- 
cal scheme i s the Total Variation Diminishing Lax-Friedrich 
scheme (see lT6tblll996T ) with a third order accurate Koren 
limiter (|Korenlll993F on the primitive variables. We use a 
uniform cylindrical grid with r/ro £ [0.2,2.4], and the full 
azimuthal direction tp £ [0, 2ir], Numerical convergence has 
been tested and two different resolutions have been used: 
(1024, 128) for most of the simulations (m = 1 to 4) but 
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Figure 2. Effective potential for a density profile given by Eg. 1131 
with x = 0-2 for m = 2 (upper panel) and m = 5 (lower panel). 
Waves can propagate only in the regions where V c g (r) < 0. 



a resolution of (1024, 256) was needed for the highest az- 
imuthal modes number m = 5 and the simulations with 
white noise perturbations. There is a null radial velocity at 
the inner and outer boundary of the grid. 

This numerical co nfiguration is very similar to the one 
of iMeheut etaH (|2012h which has proven to accurately de- 
scribe the instability, with a very good agreement with the 
analytical approach in the linear phase. For all these simu- 
lations, the square of the epicyclic frequency, k 2 , is positive 
at all radii and the disc is stable against axisymmetric per- 
turbations. 



3.2 Results 

For each of the 20 simulations, the linear growth rate of 
the instability is obtained by fitting the amplitude of the 
vorticity in the Rossby waves region as a function of time 
(see Fig. [3j . We have previously shown that this method 
gives growth rate very similar to the one ob tained by solving 
the linear equations (jMeheut et al.ll2012T j. The amplitude 
of the waves initially grows exponentially until saturation. 
As explained in section l2~2l we estimate saturation to arise 
when the waves reach a vorticity of the order of 27. This 
amplitude is plotted as a dashed line in Fig. [3] 

For each simulation, the ratio of the maximum of vortic- 
ity to the expected value, uj™ ax /2~/, is given. The saturation 
amplitude reached in the non-linear numerical simulations 
for m = 4 and m — 5. Similar to the non- linear evolution 
of the Landau damping, we also obtain oscillations in the 
waves amplitude. 

For low mode numbers, the maximum vorticity is 
slightly below the estimation. This is due to the shape of 
the vortices. Indeed the vortices with low m have elongated 
shapes as can be seen in Fig.|4]where vortices streamlines in 
the (r, tp) plane are plotted. These streamlines are computed 
with the perturbed velocity during the linear stage of the 
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Figure 3. Amplitude of the Rossby waves (vorticity) on logarithmic scale as a function of time (in units of f2 Q ). A fit of the exponential 
growth (solid line) gives the growth rate 7. The dashed line corresponds to the saturation amplitude estimated by the model. 



instability. To allow for a comparison between different m 
cases, the time slices for these plots are chosen such that the 
perturbations have about the same amplitudes. This shape 
was to be expected as the width of the Rossby wave prop- 
agation region is fixed by the width of the initial density 
bump and the length of the vortices is directly related to 
the azimuthal mode number m. Assuming a doubled cir- 
culation time for the elongated m — 2 vortices gives the 
correct saturation amplitude, as one can see in Fig. [3] for 
the [m = 2, x — 0.30] simulation. 

On the other hand, the vorticity maximum is overes- 
timated for the highest azimuthal mode number when the 
growth rate is low. See for instance the [m = 5, % = 0.15] 
simulation. We have checked that this is not related to nu- 
merical dissipation by doubling the resolution and obtain- 
ing no modification of Ld™ ax . Indeed the energy loss respon- 
sible for the low saturation amplitude may be due to the 



density waves propagating outside the Lindblad resonances 
that were not considered in the local mechanism proposed 
for growth and saturation in section [2~2l The amplitude of 
the density waves and the energy loss is the highest when 
the Lindblad resonances are close to the corotation radius 
and the width of the evanescent wave region is small. This 
regions correspond to a positive effective potential V e s in 
Fig. [21 Since the Lindblad resonances are closer to corota- 
tion for higher azimuthal mode number, the evanescent re- 
gion is narrower. Energy transmission through this region is 
then carried away by density waves. These density waves can 
then become non-linear when they reach high amplitudes. 
This may explain the lower amplitude at saturation. More- 
over these density waves are also responsible for the angular 
momentum transfer through the disc and as a result for the 
evolution of the radial structure of the disc from which the 
linear growth is computed. The distance between the two 
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Figure 4. Streamlines showing the shape of the vortices for the 
four azimuthal mode numbers considered here. For the lowest 
mode number the vortices are elongated, for the highest ones the 
vortices are circular. 



propagation region appears clearly on Fig. [5]where both the 
vorticity waves in the region of corotation and the position 
of the Lindblad resonances are visible. One can also identify 
the spiral density waves propagating inwards and outwards 
from the Lindblad resonances in the m = 5 plot. To fur- 
ther check the importance of the Lindblad resonances, we 
have performed numerical simulations with a shallower den- 
sity bump (er/rn = 0.4) but with the same azimuthal mode 
number m — 5. The density bump being shallower, the Lin- 
blad resonances are situated further away from the vortices 
and less energy is lost in the spirals waves. Fig.[5]shows that 
in this case the saturation amplitude is correctly estimated. 

The amplitude of saturation in the simulations with an 
initial white noise perturbation is more complex. Multiple 
modes grow simultaneously in the linear phase with different 
growth rates, as can be seen on Fig(7] However we still ob- 
tain an exponential growth of the total perturbation as the 
growth rates of the different modes are similar. The growth 
rate of the total perturbation (Q.21fio) is slightly smaller 
than the highest growth rate of the different modes (0.22fio 
for m = 4) due to this coexistence of multiple modes with 



m = 2 m = 5 




0.6 1.0 1.4 0.6 1.0 1.4 

r/r„ r/r 

Figure 5. Vorticity (Oo) of the flow near corotation for m2x25 
and m5x25. The white lines show the position of the Lindblad 
resonances. The width of the evanescent region between the vor- 
ticity wave and the Lindblad resonance is larger for lower m. 
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Figure 6. Amplitude of the Rossby waves (vorticity) on logarith- 
mic scale as a function of time (in units of Qq ) for the simulation 
with cr/ro = 0.4 and m = 5. A fit of the exponential growth (solid 
line) gives the growth rate 7. The dashed line corresponds to the 
saturation amplitude estimated by the model. 
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Figure 7. Amplitudes of the Rossby waves with different m's as 
a function of time. The simulation starts out with a random white 
noise, with the parameters \ = 0.3, cr/ro = 0.05, c s / (roQo) = 0.1. 
The different m component is obtained by a Fourier analysis in 
the azimuthal direction. The grow rates of the dominant modes 
are also given. 

similar growth rate. These different modes then interact in 
the non-linear phase, but one can see that eq. (|10[) still gives 
a good approximation of the saturation amplitude of the 
instability in this more realistic case. 

3.3 Fitting formula 

In order to give an estimation of the growth rate of the 
instability based on the parameters of the density bump, 
we performed simulations with initial random perturbations 
and with % £ [0.1,0.3] and a/h € [0.45,0.6]. Fitting sepa- 
rately the dependence of the growth rate on \ an( A a/h, and 
then combining, we find that 7 is approximately given by 

^(^(^- fc2 f 3 ' < lB > 

with k\ = 0.142 and k% = 0.138. Here f2o the Keplerian 
frequency at the location of the density bump, a/h and x 
specify the width and amplitude of the bump (see eg. I13[) . 
and h = c s /Qo is the disc scale height. Note that this fit- 
ting formula should be applied only in the range of param- 
eters specified above and require that the disc be stable 
against axisymmetric perturbations (i.e., the Rayleigh cri- 
terium ft 2 > 0, is satisfied). Eg. 1151 also corresponds to the 
growth rate of the fastest growing mode (m ~ 4) in the linear 
theory. The amplitude of maximum vorticity at saturation 
is given by 

3.4 3D simulation 

In the above we have considered 2D discs, neglecting their 
vertical structure. Recent study of the RWI in 3D verti- 
cally stratified discs has reveale d unexpected vertic al dis- 
placements inside the vortices (|Meheut et al.l [2010) that 
could be relevant fo r the saturation mechanism. However 
iMeheut etail l|2012l ) and [Cin| (|2012h have shown that the 
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Figure 8. Amplitude of the Rossby waves (mid-plane vorticity) 
on logarithmic scale as a function of time for the 3D simulation. 
A fit of the exponential growth gives the growth rate (solid line). 

growth rate is only slightly modified (decreased) in a fully 
stratified disc. Moreover the vertical velocity in the vortices 
is of low amplitude compared to the radial and azimuthal 
perturbed velocity. As a result, we do not expect the cir- 
culation time to be significantly modified by this vertical 
circulation. To confirm the relevance of our model for 3D 
stratified discs, we perform a numerical simulation in such 
a configuration. We choose m = 4 and x = 0.25, and the 
volume density i s chosen to have a su rface density similar to 
the 2D case (see IMeheut et all |2012h . The resulting growth 
and saturation of the instability is plotted in Fig. [8] The 
predicted maximum amplitude of the instability is correct, 
meaning that the model can be extended to 3D discs. This 
can then be also app lied to the saturation phase obtained in 
IMeheut et al.l l|2012l ) where the long-term stability of the 3D 
vortices was studied. After the saturation the vortices con- 
tinue to evolve, they tend to merge as in an inverse cascade 
and eventually to decay. 



4 SUMMARY 

We have tested numerically the non-linear saturation 
mechanism for the Rossby wave instability proposed by 
lLovelace et alj (|2009l ) . This is based on the classical results 
of particle- wave interaction well known for unstable electron 
plasma waves. We have shown that non-linear saturation of 
the RWI occurs when the trapping frequency of the parti- 
cles in the Rossby vortices equals the instability growth rate. 
We have estimated the trapping frequency by the vorticity 
of the Rossby wave; this is accurate for circular vortices but 
is only a lower limit for elongated vortices. We also discussed 
the linear coupling between the Rossby waves and the den- 
sity waves. This has the effect of decreasing the amplitude 
of the vorticity wave. Finally, we gave an estimation for the 
amplitude of the vortices at saturation based on the charac- 
teristics of the density bump. 

Our paper focused on the non-linear saturation of the 
linear instability in an initially steady disc. We did not con- 
sidered the mechanism for the formation of the pressure 
bump causing the instability. In a disc where the bump is 
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continuously sustained by some accretion processes, the disc 
is not steady and a future step is to take into account the 
evolution of the shape of the bump in the study of the in- 
stability. 
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